%%%%%%%%% Programming Language: Matlab 2019b       %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% (C) Andrew B. Abel and Stavros Panageas 2021 %%%%%%%%%%%%%%%%%%
%%%%%%%%% This code produces the first tow columns of table 4 for the paper
%%%%%%%%% "An Analytic Framework for Interpreting %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Investment Regressions in the Presence of Financial Constraints

clear all;
clc;
global rho;
global mu_H;
global mu_L;
global phi_L;
global phi_H;
global theta;
global alpha;


dt=0.01;
rho=0.07;%0.07
mu_H=0.3;%0.3
mu_L=0.1;
phi_L=-0.15;%-0.1
phi_H=0.2; %0.2
alpha=0.2;
theta=20;
tic
rng(2525)


% Preliminary steps: Initial guesses for x^ast and m
eta_1=(rho+mu_H)/phi_L;
eta_2=(rho+mu_L)/phi_H;
A=[eta_2 -mu_L/phi_H;-mu_H/phi_L eta_1];
om=eig(A);
om_1=om(2);
om_2=om(1);
x_m=max(1/(om_2-om_1)*log(om_1^2/om_2^2*(eta_2-om_2)/(eta_2-om_1)),0);
c_1=1/(om_2-om_1)*(om_2/om_1);
c_2=1/(om_1-om_2)*(om_1/om_2);
mmf=@(y) c_1*(rho+mu_L)/mu_L*(1-om_1/eta_2)*exp(om_1*(-y))+c_2*(rho+mu_L)/mu_L*(1-om_2/eta_2)*exp(om_2*(-y)) - alpha; % NOTE THAT \gamma=0
x_m_s=fzero(mmf,x_m);



value_H=@(x) c_1*exp(om_1*(x-x_m_s))+c_2*exp(om_2*(x-x_m_s));
derivative_H=@(x) om_1*c_1*exp(om_1*(x-x_m_s))+c_2*om_2*exp(om_2*(x-x_m_s));
value_L=@(x) c_1*(rho+mu_L)/mu_L*(1-om_1/eta_2)*exp(om_1*(x-x_m_s))+c_2*(rho+mu_L)/mu_L*(1-om_2/eta_2)*exp(om_2*(x-x_m_s));
derivative_L=@(x) c_1*om_1*(rho+mu_L)/mu_L*(1-om_1/eta_2)*exp(om_1*(x-x_m_s))+c_2*om_2*(rho+mu_L)/mu_L*(1-om_2/eta_2)*exp(om_2*(x-x_m_s));

disp('Solving for the value function....')
m_cand=value_L(x_m_s)/value_H(x_m_s);
sol=fminsearch(@detm, [m_cand x_m_s]); % Solve for m*, x*
sol2=detm2(sol); % determine the value functions and their derivatives

v_H=sol2(:,5);
v_L=sol2(:,3);
v_H_p=sol2(:,6); % Derivative of v_H
v_L_p=sol2(:,4); % Derivative of v_L
xt=sol2(:,1);
x_m=sol(2);

regime(1)=1;
nr_sim=1;
x(1)=0;
i=1;
disp('Simulating the model. Please wait.')
while nr_sim<100
  flag=0;
  
  
  
    if regime(i)==1
        
        value(i)=interp1(xt,v_H,x(i),'nearest','extrap');
        investment(i)=1/theta*(interp1(xt,v_H./v_H_p,x(i),'nearest','extrap')-(1+x(i)));
        
    else
        value(i)=interp1(xt,v_L,x(i),'nearest','extrap');
        investment(i)=1/theta*(interp1(xt,v_L./v_L_p,x(i),'nearest','extrap')-(1+x(i)));
        
    end
    avq(i)=value(i)-x(i);
    
    
    if regime(i)==1
        cf(i)=phi_H-investment(i)-0.5*theta*investment(i)^2;
        x(i+1)= min(x_m, x(i)+(cf(i)-(x(i)*investment(i)))*dt);
    else
        cf(i)=phi_L-investment(i)-0.5*theta*investment(i)^2;
        x(i+1)= max(0,   x(i)+(cf(i)-(x(i)*investment(i)))*dt);
    end
    
    
    
    if regime(i)==0
        if x(i+1)<=0 
            nr_sim=nr_sim+1;
            regime(i+1)=1;
            flag=1;
            x(i+1)=0;
        end
    end
    
    
    
    
    
    if regime(i)==1 
        if rand<1-exp(-mu_L*dt)
            regime(i+1)=0;
        else
            regime(i+1)=1;
        end
    else
        if flag==0
        if rand<1-exp(-mu_H*dt)
            regime(i+1)=1;
        else
            regime(i+1)=0;
        end
        end
    end
    i=i+1;
end
    
Y=investment';
X1=[avq' cf' ones(i-1,1)];
beta1=inv(X1'*X1)*X1'*Y;

disp('This is the output for Table 4 (Last two columns)')
disp('Coefficient on Average q')
beta1(1)
disp('Coefficient on Cash Flow')
beta1(2)



toc

function dydx = detm(inp)
global rho;
global mu_H;
global mu_L;
global phi_L;
global phi_H;
global theta;
global alpha;

m=inp(1);
x_star_temp=inp(2);

B_H=@(x,y) -(1/theta*(1+x) -mu_L*(y-(mu_L+rho)/mu_L));
B_L=@(x,y) -(1/theta*(1+x) -mu_H*(1/y-(mu_H+rho)/mu_H));
C_H=@(x) phi_H+0.5/theta*(1+x)^2;
C_L=@(x) phi_L+0.5/theta*(1+x)^2;
A=0.5/theta;

z_H=@(x,y) (-B_H(x,y)-sqrt(B_H(x,y).^2 -4*A*C_H(x)))/(2*A);
z_L=@(x,y) (-B_L(x,y)+sqrt(B_L(x,y).^2 -4*A*C_L(x)))/(2*A);

dydx_temp_1 = z_L(x_star_temp,m)./z_H(x_star_temp,m)*(mu_L+rho)/mu_L-m;

[t,y]=ode45(@(t,y) -y*(1/z_L(x_star_temp-t,y)-1/z_H(x_star_temp-t,y)), 0:0.0001:x_star_temp, m);

y=flip(y);
t=x_star_temp-flip(t);
vals=1/z_L(0,y(1))*(t(2)-t(1));

for j=2:length(t)
    vals=vals+1/z_L(t(j),y(j))*(t(j)-t(j-1));
end
dydx_temp_2 = alpha*exp(vals)-z_H(x_star_temp,m)*m;

dydx=dydx_temp_1^2+dydx_temp_2^2;
end


function dydx = detm2(inp)
global rho;
global mu_H;
global mu_L;
global phi_L;
global phi_H;
global theta;
global alpha;

m=inp(1);
x_star_temp=inp(2);

B_H=@(x,y) -(1/theta*(1+x) -mu_L*(y-(mu_L+rho)/mu_L));
B_L=@(x,y) -(1/theta*(1+x) -mu_H*(1/y-(mu_H+rho)/mu_H));
C_H=@(x) phi_H+0.5/theta*(1+x)^2;
C_L=@(x) phi_L+0.5/theta*(1+x)^2;
A=0.5/theta;

z_H=@(x,y) (-B_H(x,y)-sqrt(B_H(x,y).^2 -4*A*C_H(x)))/(2*A);
z_L=@(x,y) (-B_L(x,y)+sqrt(B_L(x,y).^2 -4*A*C_L(x)))/(2*A);

[t,y]=ode45(@(t,y) -y*(1/z_L(x_star_temp-t,y)-1/z_H(x_star_temp-t,y)), 0:0.0001:x_star_temp, m);

y=flip(y);
t=x_star_temp-flip(t);
vals=zeros(2,length(t));
vals(1,1)=1/z_L(0,y(1))*(t(2)-t(1));
vals(2,1)=1/z_L(0,y(1));

for j=2:length(t)
    vals(1,j)=vals(1,j-1)+1/z_L(t(j),y(j))*(t(j)-t(j-1));
    vals(2,j)=1/z_L(t(j),y(j));
end

vals2=zeros(2,length(t));
vals2(2,1)=1/z_H(0,y(1));
for j=2:length(t)
    vals2(2,j)=1/z_H(t(j),y(j));
end

v_L=alpha*exp(vals(1,:));

v_H=v_L./y';
v_L_p=v_L.*vals(2,:);
v_H_p=v_H.*vals2(2,:);

dydx=[t y v_L' v_L_p' v_H' v_H_p'];
end

